% Multiple band processing 
% takes on band after another
% bit fustrtating but just works

% latest date 28.03.2003

% code to compute MT trfn from raw ts files 
% Support RTS files created by GMS05 
% usage -> mtprocn 
% uses a suit of sub-functions.

% references-----------
% 1) Sims, 1971 (ls impedance) 2) Dekker, (predicted coherency) 1984 
% 3) Gamble et al
% 4) 1979 (rr)	4)ProcMT manual (rts file structue, biv error)
% 5) Manoj & Nagarajan, 2003 (ANN selection)
% 6) Chave et al 1989 (Jack Knife) 7) Ritter et al 1998 (Robust processing)

fllag=1;
pp=0;

while fllag==1,
   
[fid,er,file] = getbfile('\\Ibm\e\manoj\data\sgt\*.rts');

if fid == -1 | er < 1,
   fllag=0;
end;
if fllag==1,
   

disp('Rts-in -> reading time series');
[TS,fid] = rts_in(fid); % reading time series file (rts format - metronix)
fclose(fid);

ProcDef = SetDefault(TS.head.band); % setting defaults. ProcDef is a structured array
ProcDef.file1=file;

if ProcDef.ref == 1, % rr processing
   disp('Show me the reference time series !');
   [fid,er,file] = getbfile('*.rts');
if fid == -1 | er < 1,
 	return;
end;
disp('Rts-in -> reading reference time series');
  	[TSr,fid] = rts_in(fid);
  	fclose(fid);
  	disp('RrCorrectTS -> Checking overlap, trimming offshoots.. ');
  	[TS,er,len] = RrCorrectTS(TS,TSr,ProcDef); % check for overlap, trim extra lengths, pack the reference field to the local filed
  
  	if er<1,
	   fprintf('Error file mismatch\n');
	   er=-1;
   return;
   else,
   	ProcDef.file2=file;
   	clear TSr;
   end;
end;

ProcDef.len = floor(length(TS.matrix)/ProcDef.block);

ProcDef.SelStacks = ones([1 ProcDef.len]);  

if ProcDef.nband == 4 & ProcDef.DoANN == 1,
	ProcDef = SelANN(TS,ProcDef);
end;                                               
                                               
                                            
SPMatrix = GetSpectra(TS,ProcDef);
temp = size(SPMatrix); % This is required, as number of overlapped ts sections may not be just 2*len-1
ProcDef.len = temp(2);

ProcDef.Weights = ones([ProcDef.nfrq ProcDef.len]); % Weights are modified in robust averaging routine(s)

disp('AverageSPM (Non robust)-> Averaging  spectra sets....');

%SPMatrix = ScaleSPM(SPMatrix,ProcDef);
%SPMatrix = NormSPM(SPMatrix,ProcDef); % Normalize the SPMs wrt power in horz. mag fields

SPM  = AverageSPM(SPMatrix,ProcDef);

for p = 1:ProcDef.nfrq,
   SPMa.spectra(pp+p) = SPM.spectra(p);
end;
pp = pp+ProcDef.nfrq;
end;
clear SPM SPMatrix
end;

if pp>0,

	for i = 1:pp,
   	fq(i)=SPMa.spectra(i).freq;
	end;

	ProcDef.TLFreq1=fq;
	ProcDef.nfrq=length(fq);

	tf1 = transfn(SPMa,ProcDef);
   [Rhxy,Rhyx,Phxy,Phyx] = dispresph(tf1,ProcDef);
   
end;






 



